function outdir = ratheartT2fit( studydir, file, varargin )
% ratheartT2fit: perform T2 fitting for studydir.
%
% Usage: 
%   outdir = ratheartT2fit( studydir )
%   outdir = ratheartT2fit( studydir, force )
%
% Where:
%   studydir = the T2 map in the dicom folder, e.g. <pathtoratdata>\20110817\1097\T2
%   force = true/false (default false), if set to true, force redrawing
%   ROIs.
%   outdir = the directory with the output results.
% 
% Results are written to the outdir, a directory named <studydir>postproc.
%
% Based on Henk's T2Execute.m.
% Original documentation:
%
% T2mapCramerRao: Estimates parameters and errorbars for T2 relaxation 
% model = A*exp(-TE/T2)% 
%
% Two ROIs need to be drawn: First a ROI with only noise and second the ROI
% where the mapping has to be performed
%
% Output are dicom images of the Amap, R2map, R2errormap, T1map and T1errormap, saved
% in the directory of the first image of the serie. IMPORTANT: The values of the R2 and
% R2 errormap are multiplied with 10000 to get into sensible gray values.
% The errors in the errormap are the lower bounds on the standard deviations calculated with
% the Cramer-Rao lower bounds.
%
% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot, University of Antwerp, 13-8-2007.

%clear all

force = false;
if nargin > 2
    force = varargin{1}
end

roibased = 0; % 0: pixel based mapping. 1: roibased mapping.
niftioutput = 0; % 0: only .dcm output. 1: .dcm and .nii output

% size of figure for manual segmentation:
figuresize = [100 100 1000 800];

%The directory of the MR images
%studydir='D:\data\proj\ratheart\repro\20110817\1097\T2';
% now use argument

if ~exist(studydir)
    error( 'ratheartT2fit:invalidstudydir', [ studydir ' does not exist' ] );
end

% remove slashes at the end
while ( studydir(end) == '/' || studydir(end) == '\' )
    studydir = studydir(1:(end-1));
end

outdir = [ studydir 'postproc\' ]
if ~exist(outdir, 'dir')
  disp( ['making' outdir]);     
  [success, message, messageid] = mkdir(outdir);
  if ~success    
     error( 'ratheartT2fit:invalidoutputdir', message );     
  end 
end

% file='0001';

%maskdir='C:\Users\hsmit\Desktop\Data\429\026_BBFSE_TE18';
%maskfile='T2starmap0001';
%maskinfo = dicominfo([maskdir,'\',maskfile]);
%maskimage = double(dicomread(maskinfo));

numberofROIs=1;

cd(studydir)
d=dir;
nfol=length(d);
fitdata(1).folder=char(d(1+2).name);

%get the dicom info
for i=1:nfol-2
    fitdata(i).folder=char(d(i+2).name);
    info = dicominfo([studydir,'\',fitdata(i).folder,'\',file]);
    fitdata(i).image = double(dicomread(info));
    fitdata(i).TE=info.EchoTime;
    fitdata(i).TR=info.RepetitionTime;
    fitdata(i).Angle=info.FlipAngle;
    fitdata(i).info=info;
end

TE=zeros(nfol-2,1);
for k=1:nfol-2
    TE(k)=double(fitdata(k).TE);
end
% find which one has the lowest TE
[lowestTE, firstTEindex] = min(TE);
firstTEimage = fitdata(firstTEindex).image;

%determine the sigma of the noise
noisemaskfile = [outdir 'noisemask.mat'];
if ( ~exist( noisemaskfile, 'file' ) || force )
  figure( 'position', figuresize );
  imshow(fitdata(firstTEindex).image,[]);
%   title( 'draw region in background, to estimate noise level','FontSize',20,'FontWeight','bold' );
  text(35,120,'Select a region in the background','Color','yellow')
  drawnow;
  M=roipoly;
  close;

  % save this mask
  save(noisemaskfile, 'M')
else
  % load previously used noise mask file.
  disp( 'using previously generated noisemask.mat file' )
  load(noisemaskfile, 'M');  
end

for i=1:nfol-2
    D=fitdata(i).image.*M;
    sdyd(i)=std(nonzeros(D));
    maxsd=max(sdyd);    
end


clear D;

%conversion factor from rayleigh distribution
CRLBsigma=1.527*maxsd; 

%read in data to fitdata struct
cd(studydir)
d=dir;
nfol=length(d);

%when ROI is to be drawn here:
heartroifile = [outdir 'heartroi.mat'];
if ( ~exist( heartroifile, 'file' ) || force )

  figure( 'position', figuresize );  
  imshow(fitdata(firstTEindex).image,[]);    
  title( 'draw outer boundary of myocardium (end with double-click)' );
  drawnow;
  M=roipoly;
  title( 'draw inner boundary of myocardium (end with double-click)' );
  N=roipoly; % outer and inner circle of myocardium
  M = M - N;
  close;
  save( heartroifile, 'M');
else
  disp( 'using previously generated heartroi.mat file' )
  load( heartroifile, 'M');          
end
%when ROI is imported
%   M = maskimage ./ maskimage;
%do the whole image   
%   M=ones(256,256);

for i=1:nfol-2
    fitdata(i).image=fitdata(i).image.*M;
end

%initialize&allocate
TEprime=TE;
yd=zeros(size(TE,1),size(TE,2));

if ~roibased
    [row,column]=find(fitdata(2).image);
    nrvoxels=size(row,1);
end;


dims=size(fitdata(1).image);
R2=zeros(dims(1),dims(2));

dims=size(fitdata(1).image);
T2=zeros(dims(1),dims(2));
R2=T2;
A=T2;
STDCRT2=T2;
STDCRR2=T2;
STDCRA=T2;

% CRT1=T1;
% CRR1=T1;
% CRA=T1;
% CRB=T1;
% T1thresh=T1;
% EigVec=zeros(3,3,nrvoxels);
% EigVal=EigVec;

%setting fit options
opt = optimset('fminunc');
opt = optimset(opt,'LargeScale','off','DerivativeCheck','off','gradObj','on','Display','off','MaxIter',40,'Hessian','off','TolFun',1e-15,'TolX',1e-7);
s = fitoptions('Method','NonLinearLeastSquares','Lower',[-10000,-1000,-1000],'Upper',[50000,10000,10000],'Startpoint',[1000,0,0],'Maxiter',30,'TolFun',10^-9,'TolX',10^-9,'Display','off');
f = fittype('a*exp(-x*b)+c','options',s);
fitrange = 1:size(TE);

 if roibased
         yd=double(fitdata(1).image);
        
         for k=2:nfol-2
                yd=[yd double(fitdata(k).image)];
         end
         
         yd=reshape(yd,size(yd,1)*size(yd,2),1);
         TE=repmat(TE,1,size(yd,1)/(nfol-2));
         TE=TE';
         TE=reshape(TE,size(yd,1),1);
         [nonzerox,nonzeroy]=find(yd);
         yd=yd(nonzerox,1);
         TE=TE(nonzerox,1);
         nrvoxels=1;
 end
    
for i=1:nrvoxels
    disp(['Calculating pixel: ' num2str(i) '/' num2str(nrvoxels)]);
    if ~roibased   
        if fitdata(1).image(row(i),column(i)) > 0;
                for k=1:nfol-2
                    yd(k)=double(fitdata(k).image(row(i),column(i)));
                end

                [LS,ML]=T2ABCComputePar(yd,TE,0,CRLBsigma,opt,s,f, fitrange);
                [CR]=T2ABCCramerRao(ML,TE,CRLBsigma);

               % if size(find(CR>0,1))>3
               %     [EigVec(:,:,i),EigVal(:,:,i)]=eig(CR);  
               % else
               %     [EigVec(:,:,i)]=[0 0;0 0];
               %     [EigVal(:,:,i)]=[0 0;0 0]; 
               % end

                A(row(i),column(i)) = ML(1,1);
                R2(row(i),column(i)) = ML(2,1);
                T2(row(i),column(i)) = 1/ML(2,1);

                %CRA(row(i),column(i)) = CR(1,1);
                %CRR2(row(i),column(i)) = CR(2,2);
                %CRT2(row(i),column(i)) = CR(2,2)/(ML(2,1)^4);

                if CR(1,1)>0 && ~isnan(CR(1,1))
                    STDCRA(row(i),column(i)) = sqrt(CR(1,1));
                else
                    STDCRA(row(i),column(i)) = 10000;
                end
                if CR(2,2)>0 && ~isnan(CR(2,2))
                    STDCRR2(row(i),column(i)) = sqrt(CR(2,2));
                    STDCRT2(row(i),column(i)) = sqrt(CR(2,2))/(ML(2,1)^2);
                else
                    STDCRR2(row(i),column(i)) = 5;
                    STDCRT2(row(i),column(i)) = 10000;
                end

            else
                A(row(i),column(i))= 0;
                T2(row(i),column(i)) = 0;
                R2(row(i),column(i)) = 0;
 
%                 CRA(row(i),column(i)) =0;
%                 CRR2(row(i),column(i)) = 0;
%                 CRT2(row(i),column(i)) = 0;
                STDCRR2(row(i),column(i)) = 0;
                STDCRT2(row(i),column(i)) = 0;
        end  
    else
        [LS,ML]=T2ComputePar(yd,TE,0,CRLBsigma,roibased,opt, fitrange);
        [CR]=T2CramerRao(ML,TE,CRLBsigma);
        A = ML(1,1);
        T1 = 1/ML(2,1);
        R1 = ML(2,1);
        if CR(1,1)>0 && ~isnan(CR(1,1))
            STDCRA = sqrt(CR(1,1));
        else
            STDCRA = 10000;
        end
        if(CR(2,2)>0 && ~isnan(CR(2,2)))
            STDCRR1 = sqrt(CR(2,2));
            STDCRT1 = sqrt(CR(2,2))/(ML(2,1)^2); %error propagation dT/T = dR/R
        else
            STDCRR1 = 5;
            STDCRT1 = 10000;
        end
    end
end

%write result in dcm files. values*1000 to display the values in
%milliseconds, sqrt of CRR2 is to go to stds instead of vars.

%lx=nonzeros(EigVec(1,1,:)).*sqrt(nonzeros(EigVal(1,1,:)));
%ux=nonzeros(EigVec(1,2,:)).*sqrt(nonzeros(EigVal(2,2,:)));
%ly=nonzeros(EigVec(2,1,:)).*sqrt(nonzeros(EigVal(1,1,:)));
%uy=nonzeros(EigVec(2,2,:)).*sqrt(nonzeros(EigVal(2,2,:)));

if ~roibased
    dicomwrite(int16(10000*R2),[outdir,'\','R2map','.dcm'],fitdata(1).info);
    dicomwrite(int16(T2),[outdir,'\','T2map','.dcm'],fitdata(1).info);
    dicomwrite(int16(A),[outdir,'\','A2map','.dcm'],fitdata(1).info);
    dicomwrite(int16(10000*STDCRR2),[outdir,'\','R2errormap','.dcm'],fitdata(1).info);
    dicomwrite(int16(STDCRT2),[outdir,'\','T2errormap','.dcm'],fitdata(1).info);
    
    % also write the first echo time image, for use by registration:
    dicomwrite(int16(firstTEimage), [outdir,'\','FirstIM', '.dcm'], fitdata(firstTEindex).info );
    imshow(T2,[]);
end

if roibased
    scatter(TE,yd,'.b')
    hold on;
    te=floor(min(TE)):1:round(max(TE));
    yy=ML(1,1)*exp(-te*ML(2,1));
    plot(te,yy,'k');
    disp(['Results: value +- standard deviation:'])
    disp(['R2 = ' num2str(ML(2,1)) ' +- ' num2str(sqrt(CR(2,2)))])
    disp(['T2 = ' num2str(1/ML(2,1)) ' +- ' num2str(sqrt(CR(2,2))/(ML(2,1)^2))])
end;
    

if niftioutput
    
        T = [0 1 0 0;-1 0 0 0;0 0 -1 0;0 0 0 1]; %useplan(k,2).Tf;
        voxsp = sqrt(sum(T(:,1:end-1).^2,1));

        R2nii = make_nii(R2, voxsp);
        R2nii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
        R2nii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
        R2nii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
        R2nii.hdr.hist.sform_code = 1;
        save_nii(R2nii, [outdir,'\','R2mapnii'])

        T2nii = make_nii(T2, voxsp);
        T2nii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
        T2nii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
        T2nii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
        T2nii.hdr.hist.sform_code = 1;
        save_nii(T2nii, [outdir,'\','T2mapnii'])

        Anii = make_nii(A, voxsp);
        Anii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
        Anii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
        Anii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
        Anii.hdr.hist.sform_code = 1;
        save_nii(Anii, [outdir,'\','Amapnii'])

        R2errnii = make_nii(STDCRR2, voxsp);
        R2errnii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
        R2errnii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
        R2errnii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
        R2errnii.hdr.hist.sform_code = 1;
        save_nii(R2errnii, [outdir,'\','R2errormapnii'])

        T2errnii = make_nii(STDCRT2, voxsp);
        T2errnii.hdr.hist.srow_x = T(1,:)./[voxsp 1];
        T2errnii.hdr.hist.srow_y = T(2,:)./[voxsp 1];
        T2errnii.hdr.hist.srow_z = T(3,:)./[voxsp 1];
        T2errnii.hdr.hist.sform_code = 1;
        save_nii(T2errnii, [outdir,'\','T2errormapnii'])     
end
    

disp( [ 'Standard deviation range of background noise: ' num2str([min(sdyd) max(sdyd)], 'min %f  max %f ') ] );
if min(sdyd) < 0.5*maxsd    
    warning( 'Standard deviation of background noise changes more than 50 percent!');
end

%To show an A vs. R2 map with errorbars(either parallel with x and y axis
%or in the direction of the eigenvectors of the covariancematrix) enable one of the following lines: 

%errorbarxy(nonzeros(A)',nonzeros(R2)',sqrt(nonzeros(CRA))',(sqrt(nonzeros(CRR2)))',[],[],'w')
%errorbarxyoblique(nonzeros(A)',nonzeros(R2)',lx',ly',ux',uy','w')
%hold on
%scatter(nonzeros(A),nonzeros(R2),'.k')
 

